Estimate of the Cosmological Bispectrum from the MAXIMA-1 
Cosmic Microwave Background Map 
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We use the measurement of the cosmic microwave background taken during the MAXIMA-1 flight to 
estimate the bispectrum of cosmological perturbations. We propose an estimator for the bispectrum 
that is appropriate in the flat sky approximation, apply it to the MAXIMA-1 data and evaluate 
errors using bootstrap methods. We compare the estimated value with what would be expected if 
the sky signal were Gaussian and find that it is indeed consistent, with a \ 2 P er degree of freedom 
of approximately unity. This measurement places constraints on models of inflation. 
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PACS Numbers : 98.80.Cq, 98.70.Vc, 98.80. Hw 

Introduction: All theories of structure formation in the 
universe predict the properties of the probability distri- 
bution function (PDF) of cosmological perturbations. In 
all cases of interest, the PDF can be completely described 
in terms of its spatial n-point correlation functions, which 
are the expectation values of all possible products of the 
random held with itself at different points in space. Un- 
der the assumption of statistical isotropy and homogene- 
ity, it is normally more useful to characterize the PDF in 
terms of higher order moments of the Fourier transform 
of the field. Most readers are familiar with the 2-point 
moment, the power spectrum of fluctuations (Ce). Indeed 
current efforts in the analysis of Cosmic Microwave Back- 
ground (CMB) data have focused mainly on increasingly 
precise estimates of the angular power spectrum. The 
theoretical bias for this is clear: for Inflation induced 
perturbations, which is the current favourite model of 
structure formation, the statistics are Gaussian and all 
non-zero moments of order n > 2 can be expressed in 
terms of the Ci. 

In this letter we present the hrst estimate of the bis- 
pectrum of the CMB on degree, and sub-degree, angu- 
lar scales. The bispectrum is the cubic moment of the 
Fourier transform of the temperature field and it can be 
seen as a scale dependent decomposition of the skewness 
of the fluctuations (in much the same way as the Ci is 
a scale dependent decomposition of the variance of fluc- 
tuations). The bispectrum can be used to look for the 
presence of a non-Gaussian signal in the CMB sky. We 
use the data collected with the MAXIMA-1 experiment 
B to quantify the bispectrum of the CMB. The Gaus- 
sianity of this data set has already been analysed using 
complementary methods in [0 , including the methods of 



moments, cumulants, the Kolmogorov test, the % 2 test, 
and Minkowski functionals in eigen, real, Wiener-filtered 
and signal- whitened spaces. 

In the past few years, interest in the bispectrum has 
grown in the scientific community. Estimates of the bis- 
pectrum in the COBE data proved the statistic to be 
extremely sensitive to some non-Gaussian features in the 
data, be they cosmological or systematic j|; the qual- 
ity of galaxy surveys has made it possible to test for the 
hypothesis that the matter overdensity is a result of non- 
linear gravitational collapse of Gaussian initial conditions 
M. On the other hand a serious effort has been under- 
taken to calculate the expected bispectrum from various 
cosmological effects; secondary anisotropies (such as the 
Ostriker-Vishniac effect, lensing, Sunyaev-Zel'dovich ef- 
fect) 0, as well as primordial sources (such as non-linear 
corrections to inflationary perturbations or cosmic seeds) 
may lead to observable signatures in the bispectrum of 
the CMB [|-|. 

Let us establish some notation. We shall be work- 
ing in the small sky approximation where a map of the 
CMB can be considered approximately flat |I9]. The 
anisotropy of the CMB, AT(x), can then be expanded in 
terms of 2-dimensional Fourier modes as follows: 



AT(x) = 



d 2 k 

(27T) 2 



a(k)e 



ik-x 



(i) 



As stated above, the complete statistical properties of 
AT can be encoded in the expectation values of prod- 
ucts of the a(k). The power spectrum is defined to 
be (a(k)a(k')) = (27r) 2 C(fc)<5 2 (k + k'). On small an- 
gular scales, the correspondence between the flat sky 
power spectrum and the full sky angular power spectum 



is straightforward: Cg — C(k)\k=e- The bispectrum is 
defined to be 

(a(k 1 )a(k 2 )a(k 3 )) = (2 7 r) 2 B(k 1 ,k 2 ,k 3 ) ( 5 2 (k 1 + k 2 + k 3 ) 

(2) 

where the delta function constraint is a consequence of 
the assumption of statistical isotropy. 
Method: In this letter we take the approach adopted by 
Ferreira, Magueijo & Gorski in the analysis of the COBE 
4 year data || : we construct an estimator for the bispec- 
trum, apply it to the MAXIMA- 1 data and quantify its 
variance using Monte Carlo methods. The MAXIMA- 1 
experiment and dataset is described in detail in Ref |L|; 
as in P| we use a map with square pixels of 8' each. 

Given a map, we Fast Fourier Transform it and con- 
struct the following bispectrum estimator: 
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where Su.^i\ is a ring in Fourier space centered at k = 
and with radial coordinates k € \i% — At/2,li + Ai/2], 
^i,feA;A< are the number of modes which satisfy this 
condition and lZe{A} is the real part of A. For a given 
choice of £{ (with i = 1, 2, 3) we obtain an estimate of the 
bispectrum averaged in a bin of width Ag. We correct for 
the finite resolution of the experiment and the pixcliza- 
tion of the map by replacing the quantity a(k) (that is 
estimated directly from the map) by a(k)/[i?(k)14 7 (k)], 
where -B(k) and W(k) are the beam and pixel window 
functions, respectively (see O] for a detailed Fourier 
space description of the beam) . 

There are a number of approximations in our analy- 
sis. We do not discuss any systematic effects that may 
have come into play when generating the map; a de- 
tailed description of these effects is presented in [ |I2[ . 
The flat sky approximation in the estimate of the power 
spectrum is valid to within 1% for the MAXIMA- 1 100- 
square-degrees map. The fact that we are not consider- 
ing a full sky map leads to two further complications p3j . 
Firstly, there will be a finite correlation length in Fourier 
space between adjacent modes. In Maximum Likelihood 
Methods this is automatically taken into account when 
constructing the correlation matrix, but in our case we 
must take care in assessing how our results depend on the 
width of the bins, Ai , in which we estimate our bispec- 
trum. Secondly, the map we are working with does not 
have periodic boundary conditions, an essential under- 
lying assumption when performing a Fast (or Discrete) 
Fourier Transform. We correct for this by multiplying 
the map by a Welch window function which suppresses 
the mismatch at the border of the map thus reducing 
the leakage between neighboring scales in Fourier space. 
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FIG. 1. Estimate of the bispectrum of the MAXIMA-1 
CMB map. The error bars are evaluated using a Monte Carlo 
bootstrap method. Note that given the small number of sam- 
ples for the low £ components, there is a large uncertainty in 
the estimation of the error bars. 

Naturally we take this into account when estimating the 
bispectrum. Finally, the map to which we apply our es- 
timator will contain anisotropic instrumental noise, and 
one may be concerned that this may bias the estimate. 
However, given that the signal and noise are uncorre- 
cted, and the noise is Gaussian p2|, it will not affect 
our estimate of the bispectrum, merely its variance. 

Our goal in this letter is twofold. Firstly to obtain an 
estimate of the bispectrum from the data without mak- 
ing any assumptions about the statistics of the signal and 
secondly to assess how compatible our estimate of the 
bispectrum is with the assumption that the MAXIMA- 
1 data set is Gaussian. To obtain a model independent 
estimate of the bispectrum and its variance we use boot- 
strap methods H). Bootstrap methods are widely used 
in situations where one wishes to extract the statistical 
properties of a given estimator without making any as- 
sumptions about the distribution from which a sample 
has been drawn. 

One can redefine the estimator in Eq. (ph in the follow- 
ing way: divide the ring in Fourier space into six equally 
sized angular segments of width 2n/6; subdivide each 
of these segments into M — 2nl/(QAi) angular slices. 
Within each of these slices apply the estimator in Eq. (H), 
replacing Sm &£) by the corresponding set of points within 
the slice. Note that this is only applicable to the diago- 
nal components of the bispectrum Bgu (the inclusion of 
non-diagonal components will introduce correlations be- 
tween samples which will bias bootstrap estimates). In 
this way we find M approximately independent estimates 
of the But; note that A? > 2ir/ (field size) for this to be 
possible. If we find the average of these M estimates we 
recover the value one obtains by applying (||). The fact 
that we have M (almost) independent estimates puts us 
in the condition where one can apply bootstrap methods 
to estimate the distribution and consequently the vari- 



ance. We should note however that there are two limi- 
tations to this approximation. On the one hand the sky 
signal is not uniformly distributed in Fourier space, i.e. 
there may be weak correlations between different Fourier 
modes. On the other hand the noise is anisotropic and 
correlated which means that the noise covariance matrix 
is not diagonal in Fourier space. Both of these effects may 
lead to correlations between the M approximately inde- 
pendent samples but for large enough A^ they should be 
negligible. Given that the bootstrap method is the only 
non-parametric (or model independent) method which 
one can apply in this situation, we choose to neglect these 
correlations [Oj. 

Our approach to test for the Gaussianity is to generate 
10 5 Monte Carlo realizations of the MAXIMA-1 data set, 
assuming a Gaussian signal with the power spectrum of 
the best fit model to the band powers estimated in [Q. 
Note that each of these mock data sets will have a real- 
ization of the noise which obeys the full anisotropic, non- 
diagonal correlation matrix; moreover the effect of pix- 
elization and finite beam are taken into account. We then 
compare our estimate of the real data with the Gaussian 
ensemble and quantify a goodness of fit. 
Results: We present the results we have obtained 
analysing a square patch in the center of the MAXIMA- 
1 map, with 50 2 pixels. Given the dimensions of the 
map, we consider A^ — 75; these correspond to the 
bin-widths of the estimates of the Cg in 0] and lead 
to correlations of order a few percent between adja- 
cent bins. In Fig. Ill we present the diagonal elements 
(£i = £2 — £3 — £) of the estimate of the bispectrum (see 
also Table S) . Note that all values of the bispectrum are 
of order (0.001 - 0.01)Cf /2 , and the fact that B m \ l=2 iA 
is so large is mostly due to the fact that this corresponds 
to the peak value of Cg. 

The boostrap errors are evaluated from resamplings 
with replacement of the approximately independent sam- 
ples within each ring; the errors correspond to the 68% 
confidence regions with these simulated distributions. 
We find that the average bootstrap errors, ai, s over an 
ensemble of Gaussian maps to be between 4% and 8% 
lower than the true underlying variance. This bias is due 
to the correlations between adjacent samples within each 
ring. Moreover the number of approximately indepen- 
dent samples ranges from M — 2 at £ — 148 to M = 10 
at £ — 748 and one should therefore bear in mind that, 
for low £ the variance in the estimate of the bootstrap 
errors is large. 

We have performed a number of tests to evaluate how 
robust the result is on the parameters of our estimator. 
We have taken a larger patch of the MAXIMA-1 map 
and considered maps of 50 2 pixels with different locations 
within the MAXIMA-1 map. The estimated bispectra 
vary by a few percent. Alternatively we have considered 
different bin-widths (with Ag — 60 and Ag = 90) and 
found that estimates of the bispectrum vary smoothly 
and are consistent within different binnings. The use of 
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FIG. 2. Comparison with a gaussian sky. The solid 
(dashed) lines delimit the 68% (95%) confidence region de- 
termined from a Monte Carlo simulation of a Gaussian sky; 
the MAXIMA-1 noise covariance matrix was used to simulate 
realistic, anisotropic noise and the beam and pixel window 
functions were included. 

the Welch window function turns out to be essential for 
small £; this is to be expected as it should be the values 
of Bgu for low £ which are most affected by finite size 
effects. A different choice of window function (such as 
the Bartlett window function) changes the estimate of 
Bui 1^=148 and -B«f|^=224 by an order of 15% but leaves 
the remaining values of Bui unaffected. One final test 
we have undertaken was to rotate the ring considered in 
Fourier space, this way displacing the M angular slices; 
we have found that the results vary by at most 10% in 
the lowest £ bin. 

In Fig H we plot the diagonal estimate of the 
MAXIMA-1 bispectrum compared to the 68% and 95% 
contour values if the sky was indeed Gaussian. We have 
checked that our statistic is unbiased even in the presence 
of anisotropic Gaussian noise and, as can be seen, the 
MAXIMA-1 Bg seem to be consistent with the Gaussian 
assumption. The obvious way to quantify this is to use a 
goodness of fit. For the Monte Carlo realizations of the 
Gaussian sky signal, we find that most of the histograms 
of the Bi-,i £, are well approximated by Gaussians and we 

jobs 

Bf^ where B% 
and C is the covariance matrix of the estimators evalu- 
ated from the Monte-Carlo realizations. In all we have 
115 values and we find \ 2 = 130. From 10 4 realizations 
we construct the expected distribution of this x 2 '- we nn d 
that 70% of the distribution is contained to the left of the 
measured value. Even if we remove the outlier from the 
set of bins centered at £ — 224 we still find that 52% of 
the distribution lies to the left of the measured % 2 . 
Cosmological Implications: One can roughly divide the 
two possible sources of non-gaussianity in the CMB into 
primordial and late time. The latter have been exten- 
sively studied in 8 and typically give rise to non-zero 



therefore define the standard x 2 = Ylt 1 1 e' I 1 1' (^£ b 'i 
I1I2I3) ^t 1 e 2 e 3 eie'e'(^e'J'e' —B„, r „ r „) where B f , f _ f _ = 



bispectra on very small angular scales [i > 1000). We 
do not expect to find any evidence for such signatures 
in the MAXIMA-1 map. Moreover, the observed bispec- 
trum limits point source contribution to the MAXIMA 
power spectrum as it shows no significant rise at high 
£. Primordial effects may give rise to non-Gaussianity 
on degree scales and we shall focus on a few possibilities 
now. Inflation predicts almost Gaussian fluctuations to a 
very good approximation; there is however the possibility 
that second order corrections in the evolution of the in- 
flaton field may lead to mild non-Gaussianity. Komatsu 
and Spergel || have parameterized this non-linearity in 
terms of a "non-linear coupling constant" , f^L , which 
can be related to dynamical parameters in a variety of 
models of inflation. For example, /jvl — (3e — 2n) where 
e and r\ are the slow roll parameters of single field in- 
flation; one expects from slow roll models that at most 
/jvl — 0(1). An order of magnitude estimate gives 



1.1 x 10 2 /AT, 
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and AT 2 = £(£+ 1) C e /(2ir). Using the Monte Carlo 
realizations described before it is possible to estimate 
the smallest amplitude, |/jvl|, distinguishable from the 
Gaussian hypothesis; we find that |/jvl| < 944 is indistin- 
guishable from a Gaussian signal at the 95% confidence 
level. Note that the use of lower multipoles (as measured 
by COBE) should narrow this interval. A fit to the mea- 
sured values using the Gaussian covariance matrix gives 
\f NL \~9O0. (x 2 = 122). 

More exotic possibilities can be considered, such as 
for example, global topological defects. A semi-analytic 
framework exists which allows one to calculate the sta- 
tistical effects using the O(N) non-linear cr- model. Dif- 
ferent values of N will correspond to different types of 
localized objects, with, for example, N = 2 correspond- 
ing to global strings, N = 3 monopoles and N — 4 cor- 
responding to textures (taking N to infinity we recover 
gaussianity) . Verde et al J7J] (see also j!| ) have estimated 
the bispectrum and found that 
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where a = A -1 / 2 (we should point out that this expres- 
sion was derived for large angles). In what follows we 
shall extrapolate this expression to subdegree scales. The 
current sensitivity is such that models with a ~< 2.4 are 
indistinguishable from Gaussian theories; this range of a 
corresponds to any value of N. We find the best fit a to 
be a = 2.2. Current estimates of the bispectrum do not 
therefore constrain global topological defects. 

The bispectrum analysis of the MAXIMA data indi- 
cates that the data is consistent with Gaussianity. This 
reinforces the conclusions obtained in M and validates 
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TABLE I. Measured bispectrum values and corresponding 
errors. The first column has the bandwidths, the second col- 
umn the estimate of the bispectrum, the third column has an 
estimate of its variance using bootstrap methods, the fourth 
column has an estimate of its variance assuming the signal 
is Gaussian and the fifth column has its variance just due to 
noise. Columns 2-5 are in units of (/iK) 3 . 

the assumptions that go in to the data-analysis pipeline, 
namely the assumption of Gaussianity of the sky signal 
which goes into both Maximum-Likelihood and Monte- 
Carlo estimates of the power spectra. 
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